import numpy

#def face_increment(N):
#	level,pos = 0,1

#	while 1:
#		if level==N: return
#		if pos==N: level = 1; pos = 0; continue
#		yield level,pos
#		pos += 1
#		if pos + level > N: level += 1; pos = 0

def face_increment(N):
	level,pos = 1,0

	while 1:
		if level == N: return
		pos += 1
		if level + pos < N:
			yield level,pos
			continue
		level,pos = level+1,0


def fill_edge(va,vb,N):
	step_ab = (vb - va)/N
	return [va + pos*step_ab for pos in range(1,N)]


def fill_face(va,vb,vc,N):
	step_ab = (vb - va)/N
	step_ac = (vc - va)/N

	face = []
	
	for level,pos in face_increment(N):
		vec = va + level*step_ac + pos*step_ab
		face.append(vec)

	return face
	
def generate_vertices():
	phi = 0.5*(1.0 + numpy.sqrt(5.0))
	twelve_points = {}
	twelve_points['X1'] = numpy.array([0.0,phi,1.0],numpy.float)
	twelve_points['X2'] = numpy.array([0.0,-phi,1.0],numpy.float)
	twelve_points['X3'] = numpy.array([0.0,-phi,-1.0],numpy.float)
	twelve_points['X4'] = numpy.array([0.0,phi,-1.0],numpy.float)

	twelve_points['Y1'] = numpy.array([1.0,0.0,phi],numpy.float)
	twelve_points['Y2'] = numpy.array([1.0,0.0,-phi],numpy.float)
	twelve_points['Y3'] = numpy.array([-1.0,0.0,-phi],numpy.float)
	twelve_points['Y4'] = numpy.array([-1.0,0.0,phi],numpy.float)

	twelve_points['Z1'] = numpy.array([phi,1.0,0.0],numpy.float)
	twelve_points['Z2'] = numpy.array([-phi,1.0,0.0],numpy.float)
	twelve_points['Z3'] = numpy.array([-phi,-1.0,0.0],numpy.float)
	twelve_points['Z4'] = numpy.array([phi,-1.0,0.0],numpy.float)

	for point in twelve_points:
		vec = twelve_points[point]
		vec /= numpy.sqrt((vec*vec).sum())

	return twelve_points


def assemble_edges():
	vertices = generate_vertices()

	edges = []
	edges.append([vertices['Z1'],vertices['Z4']])
	edges.append([vertices['Y1'],vertices['Z4']])
	edges.append([vertices['X2'],vertices['Z4']])
	edges.append([vertices['X3'],vertices['Z4']])
	edges.append([vertices['Y2'],vertices['Z4']])

	edges.append([vertices['Z1'],vertices['Y1']])
	edges.append([vertices['Y1'],vertices['X2']])
	edges.append([vertices['X2'],vertices['X3']])
	edges.append([vertices['X3'],vertices['Y2']])
	edges.append([vertices['Y2'],vertices['Z1']])

	edges.append([vertices['X1'],vertices['Y4']])
	edges.append([vertices['Y4'],vertices['Z3']])
	edges.append([vertices['Z3'],vertices['Y3']])
	edges.append([vertices['Y3'],vertices['X4']])
	edges.append([vertices['X4'],vertices['X1']])

	edges.append([vertices['X1'],vertices['Z2']])
	edges.append([vertices['Y4'],vertices['Z2']])
	edges.append([vertices['Z3'],vertices['Z2']])
	edges.append([vertices['Y3'],vertices['Z2']])
	edges.append([vertices['X4'],vertices['Z2']])

	edges.append([vertices['X1'],vertices['Y1']])
	edges.append([vertices['Y4'],vertices['X2']])
	edges.append([vertices['Z3'],vertices['X3']])
	edges.append([vertices['Y3'],vertices['Y2']])
	edges.append([vertices['X4'],vertices['Z1']])

	edges.append([vertices['Z1'],vertices['X1']])
	edges.append([vertices['Y1'],vertices['Y4']])
	edges.append([vertices['X2'],vertices['Z3']])
	edges.append([vertices['X3'],vertices['Y3']])
	edges.append([vertices['Y2'],vertices['X4']])

	return edges


def assemble_faces():
	vertices = generate_vertices()

	faces = []
	faces.append([vertices['Z4'],vertices['Z1'],vertices['Y1']])
	faces.append([vertices['Z4'],vertices['Y1'],vertices['X2']])
	faces.append([vertices['Z4'],vertices['X2'],vertices['X3']])
	faces.append([vertices['Z4'],vertices['X3'],vertices['Y2']])
	faces.append([vertices['Z4'],vertices['Y2'],vertices['Z1']])

	faces.append([vertices['Z1'],vertices['X1'],vertices['Y1']])
	faces.append([vertices['Y1'],vertices['Y4'],vertices['X2']])
	faces.append([vertices['X2'],vertices['Z3'],vertices['X3']])
	faces.append([vertices['X3'],vertices['Y3'],vertices['Y2']])
	faces.append([vertices['Y2'],vertices['X4'],vertices['Z1']])

	faces.append([vertices['Y1'],vertices['X1'],vertices['Y4']])
	faces.append([vertices['X2'],vertices['Y4'],vertices['Z3']])
	faces.append([vertices['X3'],vertices['Z3'],vertices['Y3']])
	faces.append([vertices['Y2'],vertices['Y3'],vertices['X4']])
	faces.append([vertices['Z1'],vertices['X4'],vertices['X1']])

	faces.append([vertices['X1'],vertices['Z2'],vertices['Y4']])
	faces.append([vertices['Y4'],vertices['Z2'],vertices['Z3']])
	faces.append([vertices['Z3'],vertices['Z2'],vertices['Y3']])
	faces.append([vertices['Y3'],vertices['Z2'],vertices['X4']])
	faces.append([vertices['X4'],vertices['Z2'],vertices['X1']])
	
	return faces


def icosahedron_points(N):
	points = generate_vertices().values()

	edges = assemble_edges()
	for edge in edges:
		V1,V2 = edge
		for V in fill_edge(V1,V2,N):
			V /= numpy.sqrt((V*V).sum())
			points.append(V)

	faces = assemble_faces()
	for face in faces:
		V1,V2,V3 = face
		for V in fill_face(V1,V2,V3,N):
			V /= numpy.sqrt((V*V).sum())
			points.append(V)

	return points


def icosahedron_rotations(N):
	points = icosahedron_points(N)

	I = numpy.matrix([[1,0,0],[0,1,0],[0,0,1]],numpy.float)

	rotations = []
	for pt in points:
		RZ = I.copy()
		zx_len = numpy.sqrt(pt[0]*pt[0] + pt[2]*pt[2])
		cos_entry = zx_len
		sin_entry = pt[1]

		RZ[0,0] = cos_entry
		RZ[0,1] = -sin_entry
		RZ[1,0] = sin_entry
		RZ[1,1] = cos_entry


		RY = I.copy()

		if zx_len > 0.0:
			cos_entry = pt[2]/zx_len
			sin_entry = pt[0]/zx_len

			RY[2,2] = cos_entry
			RY[2,0] = -sin_entry
			RY[0,2] = sin_entry
			RY[0,0] = cos_entry

		rotations.append(RY*RZ)

	return rotations
